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Methods  to  Determine  Wideband  System  Response  Utilizing 
only  Narrowband  Information 


1.0  Introduction 


In  a  host  of  applications  in  engineering,  it  is  necessary  to  obtain  information  about 
a  system  over  a  broad  range.  In  most  cases  it  is  not  possible  to  evaluate  the  parameter 
of  interest  in  closed  form.  However,  either  theoretical  or  experimental  data  is  available 
in  a  narrow  band.  Generation  of  the  data  over  the  broadband  is  not  possible  or  may 
be  extremely  time  consuming.  The  principle  of  analytic  continuation  is  utilized  to 
extrapolate/interpolate  the  data  over  a  wide  band.  The  method  of  Cauchy  [1]  has 
been  chosen  to  carry  out  the  analytic  continuation. 

The  Cauchy  method  deals  with  approximating  a  function  by  a  ratio  of  two  poly¬ 
nomials.  Given  the  value  of  the  function  and  its  derivatives  at  a  few  points,  the  order 
of  the  polynomials  and  their  coefficients  are  evaluated.  Once  the  coefficients  of  the 
two  polynomials  are  known,  they  can  be  used  to  generate  the  parameter  over  the 
entire  band  of  interest. 

In  this  chapter,  Cauchy’s  method  is  utilized  to  generate  broadband  currents  on 
a  body  from  which  its  Radar  Cross  Section  (RCS)  is  calculated.  This  is  done  from 
narrowband  calculations  of  the  currents.  Particularly  in  the  Method  of  Moments  [2] 
generation  of  the  response  at  each  frequency  point  is  very  time  consuming.  However, 
the  current  and  its  derivatives  with  respect  to  frequency  can  be  calculated  at  a  few 
points  using  the  Method  of  Moments.  Then  Cauchy’s  method  can  be  used  to  extrap¬ 
olate/interpolate  the  current  over  a  broad  frequency  range  from  which  the  RCS  can 
be  calculated. 

Another  example  is  the  characterization  of  optical  systems.  The  time  required  to 
evaluate  a  response  over  a  broad  range  of  the  size  parameter  would  be  prohibitive. 
The  Cauchy  method  can  be  used  to  generate  the  response  of  interest  over  a  broadband 
from  the  value  of  the  function  at  some  discrete  points. 

The  Cauchy  method  has  wide  application.  A  third  example  is  in  the  area  of  filter 
analysis.  In  the  laboratory  it  is  not  always  possible  to  make  accurate  broadband 
measurements.  This  problem  is  especially  severe  in  the  case  of  measuring  the  transfer 
function  of  a  filter  in  the  stop  band.  The  signal  to  noise  ratio  may  be  too  low  to  be 
confident  about  the  measurements  of  filter  characteristics.  Here  the  Cauchy  method 
could  be  used  to  generate  broadband  information  from  measured  narrowband  data. 

Yet  another  area  of  application  for  the  Cauchy  method  is  that  of  device  character¬ 
ization.  A  very  useful  tool  in  computer  aided  design  would  be  an  online  description 
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of  the  characteristics  of  many  devices.  But,  since  each  device  could  be  used  under 
different  operating  conditions,  each  with  its  own  frequency  characteristic,  the  memory 
required  to  describe  all  devices  would  be  prohibitive.  Here  the  Cauchy  method  could 
be  used  to  generate  broadband  information  while  storing  the  measured  data  at  only 
a  few  points. 

The  choice  of  polynomial  orders  is  restricted  by  the  information  we  have.  While 
it  is  true  that  the  more  information  we  have  the  higher  we  can  choose  the  orders, 
this  is  not  always  desirable.  In  filter  analysis  especially,  the  choice  of  the  order  of  the 
polynomials  is  very  important. 

In  this  work  the  Cauchy  technique  is  used  to  solve  the  above  problems.  In  each  of 
the  cases  mentioned  above  the  Cauchy  technique  would  save  a  significant  amount  of 
program  execution  time  or  computer  memory  while  still  producing  accurate  results 
over  broadband  frequencies.  The  method  is  tested  and  numerical  results  are  presented 
along  with  two  examples  of  the  method  as  a  time-saving  device. 


2.0  The  Cauchy  Method 


Consider  a  system  function  H{s).  The  objective  is  to  approximate  H{s)  by  a  ratio 
of  two  polynomials  A(s)  and  B{s)  so  that  H{s)  can  be  represented  by  fewer  variables. 


Consider 


H{s) 


A{s) 

B{s) 


Etohs'^ 


(2.1) 


Here  the  given  information  could  be  the  value  of  H{s)  and  its  Nj  derivatives  at  some 
frequency  points  Sj,  j  =  1,...  J.  If  H^{sj)  represents  the  derivative  of  H{s)  at 
point  s  =  Sj,  the  Cauchy  problem  is: 


Given  H^^^Sj)  for  n  =  0, . . .  Nj,  j  =  1, . . .  J,  find  P,  Q,  {a*;,  k  —  0,. . .  P],and{bk,  k  — 
0,...Q}. 

The  solution  for  {a*}  and  {6*}  is  unique  if  the  total  number  of  samples  is  greater 
than  or  equal  to  the  total  number  of  unknown  coefficients  P  -b  Q  +  2  [1],  i.e. 

N  =  Y^{Nj  +  l)>P  +  Q  +  2. 

i=i 

By  enforcing  the  equality  in  equation  (2.1)  one  obtains 

A{s)  =  H{s)B{s)  (2.2) 
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Differentiating  the  above  equation  n  times,  and  evaluating  the  expressions  at  point 
Sj,  results  in  the  binomial  expansion, 


t=0 

(2.3) 

where. 

nr 

(n-f)!*!’ 

n!  represents  the  factorial  of  n. 

Using  the  polynomial  expansions  for  A(5)  and  B{s),  equation  (2.3) 

can  be  rewritten 

as 

P  Q 

E  ^{hn),kak  =  E 

fc=0  k=0 

(2.4) 

where 

(2.5) 

B(i,o),.  =  E-a  -  0 

(2.6) 

n  =  0, 1,  ....Nj,  j  =  1,  ..J,  where  u{k)  =  0  for  A:  <  0  and  =  1  otherwise. 


Define, 


A  = 


(2.7) 


B  = 


•  •  •  ^(j,n),Q 


(2.8) 


T 

[a]  =  ao,  0^2)  ••  •  ]  (2.9) 

[(,)=(  to, 61,62,.. (2.10) 

The  order  of  matrix  A  is  iV  x  (P  +  1)  and  that  of  B  is  A  x  (Q  +  1). 

Then,  equation  (2.4)  becomes 

[A]a  =  [B]6  (2.11) 
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or 


(2.12) 


[a|-b1[;J=o 

For  ease  of  notation,  define  [C]  =  [A|  B].  C  is  of  order  N  x  {P  +  Q  +  2).  A 

Singular  Value  Decomposition  (SVD)  of  the  matrix  C  will  give  us  a  gauge  of  the 
required  values  of  P  and  Q  [3].  A  SVD  results  in  the  equation 

[U][S][Vp[;]=0  (2.13) 

where  the  superscript  H  indicates  the  matrix  Hermitian  i.e.  the  conjupte  transpose 
of  the  matrix.  The  matrices  U  and  V  are  unitary  matrices  and  S  is  a  diagonal  matrix 
with  the  singular  values  of  C  in  descending  order  as  its  entries.  The  columns  of  U 
are  the  left  eigenvectors  of  C  or  the  eigenvectors  of  CC^.  The  columns  of  V  are 
the  right  eigenvectors  of  C  or  the  eigenvectors  of  C^C.  The  singular  values  are  the 
square  roots  of  the  eigenvalues  of  the  matrix  C^C.  Therefore,  the  singular  values  of 
any  matrix  are  real  and  positive.  The  number  of  nonzero  singular  values  is  the  rank 
of  the  matrix  in  equation  (2.12)  and  so  gives  us  an  idea  of  the  information  in  this 
system  of  simultaneous  equations.  If  R  is  the  number  of  nonzero  singular  values,  the 
dimension  of  the  right  null  space  of  C  is  P  +  Q  +  2  -  i?.  Our  solution  vector  belongs  to 
this  null  space.  Hence  to  make  this  solution  unique,  we  need  to  make  the  dimension 
of  this  null  space  1  so  that  only  one  vector  defines  this  space.  Hence  P  and  Q  must 
satisfy  the  relation 

R  +  l  =  P  +  Q  +  2  (2.14) 

Hence,  the  solution  algorithm  must  include  a  method  to  estimate  R.  This  is  done  by 
starting  out  with  the  choices  of  P  and  Q  that  are  higher  than  can  be  expected  for  the 
system  at  hand.  Then  we  get  an  estimate  for  R  from  the  number  of  non-zero  singular 
values  of  the  matrix  C.  Now,  using  equation  (2.14)  we  get  better  estimates  for  P  and 
Q.  Letting  P  and  Q  stand  for  these  new  estimates  of  the  polynomial  orders,  we  can 
recalculate  the  matrices  A  and  B.  Therefore,  we  come  back  to  the  relation 


Many  methods  to  solve  equation  (2.15)  are  well  documented  [3].  For  reasons 
indicated  below,  we  choose  the  method  of  Total  Least  Squares  (TLS) 

The  usual  approach  to  solve  equation  (2.15)  is  that  of  Least  Squares  (LS).  In  this 
approach,  the  equation  is  rewritten  as: 

[Cf[C]  “  =0  (2.16) 

The  solution  for  ?  is  taken  as  the  eigenvector  corresponding  to  the  zero  eigen- 
b 

value  of  the  resulting  matrix.  However,  as  we  have  seen,  it  is  important  to  limit  the 
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rank  of  the  null  space  of  the  nnatrix  [C]  to  one.  But,  this  approach  has  an  extra  step 
of  a  matrix  multiplication.  Also,  since  the  eigenvalues  are  not  sorted,  it  is  additional 
work  to  find  the  number  of  non-zero  eigenvalues. 

A  better  approach  would  be  the  Total  Least  Squares(TLS)  [4].  Since  the  elements 
of  the  matrix  C  are  affected  by  measurement  errors  and  noise,  the  optimal  solution 
of  this  equation  must  take  into  account  the  effect  of  the  noise  in  the  matrix  elements. 
The  LS  approach  fails  to  do  so.  As  shown  in  [4],  the  Total  Least  Squares  is  the 
optimal  approach  when  the  matrix  is  affected  by  noise. 

[C]  is  a  rectangular  matrix  with  more  rows  than  columns.  Another  SVD  of  this 
matrix  brings  us  back  to  the  equation 


[U][S][V] 


=  0 


(2.17) 


By  the  theory  of  the  TLS,  the  solution  is  proportional  to  the  last  column  of  the 
matrix  V.  Since  any  constants  of  proportionality  cancel  out  while  dividing  the  two 
polynomials,  we  can  choose 


a 

b 


=  [V]p+Q+2 


(2.18) 


3.0  Applications  of  the  Cauchy  Method 


3.1  The  Method  of  Moments 

The  Method  of  Moments  yields  remarkably  accurate  solutions  to  integral  equations 
arising  in  electromagnetic  scattering  and  radiation  problems.  It  approximates  the 
interactions  of  complicated  bodies  with  a  set  of  smaller,  easily  solvable  interactions. 
The  currents  are  approximated  by  a  linear  combination  of  some  known  basis  functions. 
The  problem  then  of  extrapolating  the  current  density,  as  a  function  of  frequency, 
reduces  to  finding  the  coefficients  in  the  linear  combination.  This  approach  allows 
the  problem  to  be  written  as  a  matrix  equation  with  the  unknown  coefficients  as  the 
solution  to  the  equation.  The  Method  of  Moments  finds  it  greatest  advantage  in  the 
widespread  use  of  the  computer.  But,  its  major  limitation  is  that  a  large  matrix 
equation  has  to  be  solved  at  every  frequency  point  of  interest.  If  a  large  system  is  to 
be  studied,  the  program  execution  time  may  be  as  long  as  days. 

The  Cauchy  method  can  partially  solve  this  problem.  The  Method  of  Moments 
program  could  generate  information  over  a  limited  band  from  which  the  Cauchy 
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method  could  generate  broadband  information. 

3.1.1  Interfacing  with  the  Method  of  Moments 


The  Cauchy  method  can  easily  be  incorporated  as  part  of  a  Method  of  Moments 
program.  The  Method  of  Moments  converts  a  linear  operator  equation  into  a  matrix 
equation  of  the  form 

m  =  [z][i]  (3.1) 

Here,  [7]  is  the  vector  of  coefficients  in  the  representation  of  the  current  as  a  linear 
combination  of  basis  functions.  [V]  is  the  known  excitation  to  the  system,  while  [Z] 
is  the  matrix  that  describes  the  interaction  of  the  currents  and  the  excitation. 

Differentiating  the  above  equation  with  respect  to  frequency  results  in  a  binomial 
expansion. 

[1/]'  =  [Z]'[I]  +  [Z][I]' 

=*•  [/]'  =  [21"‘ 


In  general, 


[Ip  =  [zp 


[Vp  = 

1=1 

[vp  -j2^Ci[z]^"-pp 

i=l 


(3.3) 


In  the  above  equations,  is  the  vector  with  each  element  of  [V]  differentiated 

with  respect  to  frequency  n  times.  Similarly,  \Zp  is  the  matrix  generated  by  differ¬ 
entiating  each  element  of  the  matrix  [Z\  with  respect  to  frequency  n  times. 

Hence,  using  a  Method  of  Moments  program,  we  can  generate  all  the  information 
needed  to  apply  the  Cauchy  method.  The  use  of  derivative  information  saves  execu¬ 
tion  time  because  no  new  matrix  inversions  are  required  to  generate  the  additional 
information.  Hence,  evaluation  of  a  derivative  at  a  frequency  point  required  0{N  ) 
operations  as  opposed  to  solving  for  the  currents  at  a  frequency  point  which  takes 
0{N^)  operations.  Each  element  of  the  solution  current  ([/])  vector  is  treated  as  our 
function  H{s).  Given  the  current  and  its  derivatives  at  some  frequency  points,  we 
can  use  the  Cauchy  method  to  approximate  the  current  at  many  more  points. 

3.1.2  Numerical  Examples 


To  test  the  Cauchy  method.  Radar  Cross  Sections  {RCS)  of  five  different  perfectly 
conducting  three  dimensional  bodies  were  calculated  over  wide  frequency  bands.  A 
program  to  evaluate  the  currents  on  an  arbitrarily  shaped  closed  or  open  body  using 
the  Electric  Field  Integral  Equation  and  triangular  patching  as  described  in  [5]  was 
used.  The  triangular  patching  approximates  the  geometry  of  the  surface  of  the  body 
with  a  set  of  adjacent  triangles.  The  program  then  uses  these  currents  to  evaluate 
the  RCS  of  the  body.  It  was  modified  to  also  calculate  the  first  four  derivatives  of  the 
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currents  with  respect  to  frequency.  This  information  was  used  as  input  to  a  Cauchy 
subroutine.  The  original  Method  of  Moments  program  was  used  to  calculate  the  RCS 
without  the  Cauchy  method.  The  two  RCS  plots  were  compared  to  show  the  accuracy 
of  the  Cauchy  method. 

The  bodies  chosen  were  a  sphere,  a  square  plate,  a  disk,  a  concave,  and  a  convex 
hemisphere.  In  all  cases  the  currents  and  their  first  four  derivatives  were  evaluated  at 
five  frequency  points.  Hence,  the  total  information  allows  a  maximum  of  5  x  (4  +  1)  = 
25  coefficients  combined  in  the  two  polynomials  of  equation  (2.1).  In  the  application 
of  Cauchy  method  to  the  Method  of  Moments,  it  was  found  that  no  singular  values 
of  the  original  matrix  [yl|  —  B]  are  zero.  This  is  to  be  expected,  since,  the  current, 
as  a  function  of  frequency,  is  not  a  ratio  of  two  polynomials.  Hence,  the  higher 
the  polynomial  orders  we  choose,  the  more  accurate  the  approximation  would  be. 
Therefore,  in  this  application,  the  step  of  estimating  /?,  P,and  Q,  in  equation  (2.14) 
is  bypassed.  Given  the  25  samples,  the  numerator  polynomial  was  of  order  11  while 
the  denominator  was  a  polynomial  of  order  12.  Physically  we  know  that  for  the 
polynomial  approximation  to  be  stable,  the  numerator  polynomial  must  be  of  lower 
order  than  the  denominator  polynomial. 

The  motivation  to  apply  the  Cauchy  method  to  the  Method  of  Moments  is  to 
save  program  execution  time.  To  get  an  idea  of  how  much  time  can  be  saved,  the 
program  was  timed  for  two  of  the  above  bodies  and  compared  to  the  original  Method 
of  Moments  program.  The  two  bodies  chosen  were  the  sphere  and  the  plate. 

In  the  first  example  a  sphere  of  radius  0.3m  was  analyzed.  The  sphere  was  tri- 
angularized  using  182  nodes  and  540  edges.  Because  the  sphere  is  a  closed  object, 
this  results  in  540  unknowns  in  the  expansion  of  the  current  in  terms  of  the  basis 
functions.  The  currents  on  the  sphere  and  its  first  four  derivatives  with  respect  to 
frequency  were  evaluated  at  five  frequency  points.  The  points  chosen  were  in  the 
range  A  =  0.30m  and  A  =  0.84m  at  a  spacing  of  0.135m.  Using  this  information  and 
the  Cauchy  method  the  current  on  the  sphere  was  calculated  for  51  points  in  the  same 
frequency  range.  Using  these  currents  the  RCS  of  the  sphere  was  calculated  at  the  51 
frequency  points.  The  time  taken  for  this  calculation  is  compared  to  the  time  taken 
by  the  original  Method  of  Moments  program  to  evaluate  the  RCS  at  five  frequency 
points  in  the  same  range. 

Method  of  Moments  (5  points)  :47mins.  56secs. 

Cauchy  Method  (51  points)  :57mins.  57secs. 


To  generate  the  same  information  at  51  points  the  Method  of  Moments  program 
would  take  approximately  8hrs.8mins. 


In  Figure  3.1  we  see  the  results  of  applying  the  Cauchy  method  to  the  evaluation 
of  the  RCS  of  a  sphere.  Here  the  RCS  is  plotted  over  a  decade  bandwidth.  This 
bandwidth  was  broken  up  into  3  ranges: 
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0.6m  <  A  <  1.0m 
1.0m  <  A  <  1.8m 
1.8m  <  A  <  6.077Z 

In  each  of  the  three  ranges  the  current  and  its  first  four  derivatives  were  eval¬ 
uated  at  five  equally  spaced  points  using  the  Method  of  Moments  program.  Using 
this  information  the  polynomials  in  equation  (2.1)  were  formed.  This  rational  poly¬ 
nomial  was  used  to  evaluate  the  current  at  51  points  in  each  range.  Also,  the  original 
Method  of  Moments  program  was  used  to  calculate  the  currents  at  a  few  points  in 
the  decade  bandwidth.  The  currents  were  used  to  calculate  the  RCS  of  the  sphere  in 
this  bandwidth.  As  can  be  seen  from  the  figure,  the  agreement  between  the  results 
from  the  use  of  the  Cauchy  program  and  the  original  Method  of  Moments  program 
is  excellent. 

As  a  second  example,  a  square  plate  of  side  0.3m  was  analyzed.  The  plate  was 
triangularized  using  169  triangle  nodes  and  456  edges.  In  the  Method  of  Moments 
formulation  the  nodes  on  the  boundary  of  an  open  object  are  not  unknowns.  Hence, 
the  number  of  unknowns  in  this  case  was  only  408.  The  procedure  followed  is  similar 
to  the  analysis  of  the  sphere.  Here  the  five  frequency  points  chosen  were  in  the  range 
A  =  0.15m  and  A  =  0.30m  at  intervals  of  0.0375m.  Using  this  information  and  the 
Cauchy  method  the  currents  on  the  plate  were  evaluated  at  201  frequency  points.  The 
time  taken  for  this  calculation  is  compared  to  the  time  taken  by  the  original  Method 
of  Moments  program  to  evaluate  the  RCS  at  five  frequency  points  in  the  same  range. 

Method  of  Moments  (5  points)  :21mins.  50secs. 

Cauchy  Method  (201  points)  :27mins.  47secs. 


To  generate  the  same  information  at  201  points  the  Method  of  Moments  program 
would  take  approximately  14hrs.38mins. 


All  programs  were  executed  on  an  IBM  RS6000  platform  running  AIX. 

Figure  3.2  shows  the  application  of  this  technique  to  the  evaluation  of  the  RCS 
of  a  plate.  Again,  to  evaluate  the  RCS  over  a  decade  bandwidth,  three  intervals 
were  chosen  and  the  two  polynomials  of  equation  (2.1)  formed  in  each  interval.  The 
numerator  polynomial  had  order  11  while  the  denominator  polynomial  had  order  12. 
The  intervals  chosen  were: 


0.1m  <  A  <  0.15m 
0.15m  <  A  <  0.3m 
0.3m  <  A  <  1.0m 

The  rational  polynomial  was  used  to  evaluate  the  currents  at  201  points  in  each  range. 
The  original  Method  of  Moments  program  was  used  to  evaluate  the  RCS  of  the  plate 
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in  this  decade  bandwidth.  As  can  be  seen  from  the  figure,  the  agreement  between  the 
two  results  is  excellent. 

The  third  example  is  a  disk  of  radius  0.3m.  The  disk  was  triangularized  using 
142  nodes  and  460  edges.  Of  these  only  440  were  interior  nodes.  Figure  3.3  shows 
the  RCS  of  the  disk  over  a  decade  bandwidth.  Here  too,  the  decade  bandwidth  was 
broken  up  into  three  intervals  and  polynomials  of  order  11  and  12  formed  in  each 
interval.  The  rational  polynomial  was  used  to  evaluate  the  currents  and  then  the 
RCS  of  the  disk  at  51  frequency  points  in  each  range.  The  intervals  chosen  were: 


0.6m  <  A  <  2.4m 
2.4m  <  A  <  4.2m 
4.2m  <  A  <  6.0m 

Figures  3.4(a)  and  3.4(b)  show  the  results  of  the  final  example  of  the  Cauchy 
technique  applied  to  the  Method  of  Moments.  The  RCS  of  a  convex  and  a  concave 
hemisphere  was  calculated  over  a  decade  bandwidth.  Figure  3.4(a)  shows  the  RCS 
of  a  convex  hemisphere  while  Figure  3.4(b)  shows  the  RCS  of  a  concave  hemisphere. 
The  radius  of  both  hemispheres  was  0.3m.  The  convex  hemisphere  had  257  nodes 
and  736  edges.  This  resulted  in  a  problem  with  704  unknowns.  The  concave  hemi¬ 
sphere  had  316  nodes  and  910  edges.  Of  these  875  were  interior  nodes.  The  decade 
bandwidth  was  broken  into  the  following  ranges: 


0.6m  <  A  <  1.0m 
1.0m  <  A  <  2.6m 
2.6m  <  A  <  6.0m 

As  in  the  case  of  the  disk,  for  the  hemispheres  too,  the  two  polynomials  of  equation 
(2.1)  were  formed  in  each  of  the  three  ranges.  In  both  cases  the  Method  of  Moments 
program  evaluated  the  currents  and  its  first  four  derivatives  with  respect  to  frequency 
at  five  points  in  each  range.  This  information  was  used  by  the  Cauchy  subroutine 
to  approximate  the  currents  at  51  points  in  each  range  from  which  the  RCS  of  the 
hemispheres  were  calculated  at  51  points  in  each  range.  Also,  the  original  Method 
of  Moments  program  was  used  to  calculate  the  RCS  over  the  decade  bandwidth.  As 
can  be  seen  from  Figure  3.4,  the  agreement  in  each  case  is  excellent. 

3.2  Optical  Computations 

The  calculation  of  either  the  scattering  efficiency  or  the  intensity  is  highly  compu¬ 
tationally  intensive.  If  these  parameters  are  desired  over  a  broad  range  and  at  finely 
spaced  points  of  the  size  parameter,  the  time  required  for  the  calculations  could  be 
prohibitive.  The  Cauchy  method  would  solve  this  problem  by  needing  the  calculations 
to  be  done  at  a  much  coarser  spacing  and  interpolate  the  parameter  of  interest. 
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(a)  Convex  Hemisphere. 


Figure  3.4:  Radar  Cross  Section  of  a  Hemisphere 


This  application  was  tested  on  the  scattering  efficiency  of  a  sphere  as  a  function 
of  size  parameter[6].  The  sphere  had  an  index  of  refraction  of  2.0.  The  original  data 
calculated  the  scattering  efficiency  at  a  spacing  0.002  in  the  size  parameter.  The  range 
of  the  size  parameter  was  from  7.0  to  8.0.  Hence,  the  original  data  had  501  points. 
The  Cauchy  program  needed  a  spacing  of  0.01  in  the  size  parameter  to  accurately 
calculate  the  scattering  efficiency  of  the  sphere  at  the  original  501  points.  This  cuts 
down  the  program  execution  time  by  a  factor  of  5.  The  input  to  the  Cauchy  program 
is  shown  in  Figure  3.5(a). 

Because  all  computer  calculations  suffer  from  round-off  error,  most  of  the  singular 
values  returned  from  the  SVD  subroutine  are  not  exactly  zero.  The  choice  of  the 
threshold  was  such  that  a  singular  value  was  considered  zero  if  it  was  18  orders  of 
magnitude  lower  than  the  largest  singular  value.  This  choice  of  the  cutoff  was  chosen 
by  the  limited  precision  of  the  computer  arithmatic  i.e.  the  limited  dynamic  range  of 
the  computer.  Using  this  threshold,  only  72  singular  values  could  be  considered  not 
equal  to  zero.  Hence,  the  estimate  for  R  is  72.  Using  this  and  equation  (2.14),  the 
choice  for  the  polynomial  orders  was  reduced  to  35  for  the  numerator  and  36  for  the 
denominator.  Using  these  polynomials  the  scattering  efficiency  was  calculated  at  the 
original  501  points. 

Figure  3.5(b)  shows  the  results  of  the  application  of  the  Cauchy  method  to  optical 
computations.  The  dotted  line  represents  the  original  data  while  the  unbroken  line  the 
interpolated  data.  As  can  be  seen,  the  two  plots  are  nearly  visually  indistinguishable. 
Also,  even  though  the  input  data  to  the  program  did  not  have  of  the  peaks  of  the 
Scattering  Parameter,  the  Cauchy  program  was  able  to  reproduce  them. 

3.3  Filter  Analysis 

The  Cauchy  method  can  also  be  used  in  analysis  of  filters  over  broad  frequency 
ranges.  A  filter  response  is  a  ratio  of  two  polynomials  and  hence  lends  itself  easily 
to  the  use  of  a  Cauchy  program.  This  has  practical  application  to  the  problem  of 
generating  the  stop  band  response  given  the  pass  band  response  or  the  reverse  i.e. 
generating  the  pass  band  response  given  some  data  from  the  stop  band. 

A  filter  transfer  function  (521)  was  measured  using  a  network  analyzer  at  frequency 
points  in  and  out  of  the  filter  passband.  The  filter  had  its  3dB  points  at  4.98GHz 
and  6.61GHz.  Hence,  the  filter  had  a  passband  of  1.63GHz  with  a  center  frequency 
of  5.80GHz.  The  filter  response  was  measured  at  415  equally  spaced  points  in  the 
frequency  range  4.3lGHz  to  7.42GHz. 

In  the  first  application,  the  response  over  the  entire  band  of  measurement  was  re¬ 
covered  using  mostly  passband  information.  51  equally  spaced  points,  in  the  frequency 
range  4.79GHz  to  6.96GHz,  were  chosen  as  input  to  a  Cauchy  program.  Because  we 
are  dealing  with  measured  data,  we  do  not  have  any  information  about  the  derivative 
of  the  transfer  function  with  respect  to  frequency. 

The  threshold  was  chosen  such  that  a  singular  value  was  considered  zero  if  it 
was  14  orders  of  magnitude  lower  than  the  largest  singular  value.  After  the  program 
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(a)  Input  to  the  Cauchy  program 


(b)  Original  and  Interpolated  Data 


Figure  3.5:  Scattering  efficiency  as  a  function  of  size  parameter 
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checked  for  the  number  of  non-zero  singular  values,  the  estimate  for  R  was  16.  The 
order  of  the  numerator  polynomial  was  chosen  to  be  7  while  that  of  the  denominator 
polynomial  was  chosen  to  be  8. 

In  Figures  3.6(a)  and  3.6(b)  the  results  from  the  Cauchy  program  are  described. 
Figure  3.6(a)  shows  the  magnitude  response  while  Figure  3.6(b)  shows  the  phase 
response  of  the  filter.  As  is  often  the  case  with  filters,  the  magnitude  response  was 
considered  more  important.  Hence,  the  phase  response  was  allowed  to  show  a  poor 
agreement  so  as  to  maximize  the  agreement  of  the  magnitude  response.  If  a  10% 
error  in  the  magnitude  were  acceptable,  the  extrapolation  is  valid  for  0.39GHz.  This  is 
6.7%  of  the  center  frequency  and  23.9%  of  the  bandwidth.  For  frequencies  beyond  the 
passband,  the  extrapolation  is  accurate  within  10%  up  to  7.42GHz,  the  frequency  till 
which  data  was  available.  Hence,  we  have  generated  accurate  data  over  a  bandwidth 
of  3.32GHz  starting  with  data  over  a  bandwidth  of  2.16GHz. 

In  the  second  application,  data  from  the  stop  band  and  a  little  from  the  passband 
was  used  to  interpolate  into  the  passband.  Here  too  the  choice  of  threshold  is  very 
important.  Using  the  same  threshold  as  in  the  first  application,  the  estimate  of  R 
remained  the  same.  Hence,  in  this  case  too,  the  numerator  polynomial  had  order 
7  while  the  denominator  had  order  8.  In  this  case,  23  equally  spaced  points  from 
4.31GHz  up  to  5.35GHz  and  28  equally  spaced  points  from  6.20GHz  to  7.42GHz  were 
used  to  interpolate  into  the  passband.  This  represents  on  interpolation  of  0.85GHz, 
which  is  14.6%  of  the  center  frequency  or  52.7%  of  the  bandwidth.  Figures  3.7(a) 
and  3.7(b)  show  the  results  of  this  application.  Figure  3.7(a)  is  the  reconstructed 
magnitude  response  and  Figure  3.7(b)  is  the  reconstructed  phase  response.  Again, 
since  more  attention  is  paid  to  the  magnitude  response,  the  phase  response  shows 
poorer  agreement  with  the  true  response. 

In  both  figures  the  dotted  line  represents  the  measured  data  and  the  continuous 
line  the  results  of  the  Cauchy  program. 

3.4  Device  Characterization 

An  application  of  the  Cauchy  method  is  in  creating  a  database  of  many  devices 
working  in  varying  operating  conditions.  The  Cauchy  program  would  require  the 
value  of  a  parameter  at  a  few  frequency  points  and  use  this  information  to  evaluate 
the  paxameter  over  a  wide  frequency  band.  Over  many  devices,  and  their  operating 
conditions,  this  would  yield  significant  savings  in  memory  requirements. 

To  test  this  application  the  Y-parameters  of  a  UM  PHEMT  were  measured  over 
the  range  of  1.0-40.0GHz.  Just  five  of  these  points  were  used  as  input  to  the  Cauchy 
program.  The  points  chosen  were  at  the  frequency  points  l.OGHz,  lO.OGHz,  20.0GHz, 
30.0GHz,  and  40.0GHz.  This  resulted  in  a  numerator  polynomial  of  order  1  and 
denominator  polynomial  of  order  2.  Here  again,  the  step  of  estimating  R,  P,  and  Q 
is  bypassed.  Figure  3.8(a)  shows  the  magnitude  (IFjil)  reconstructed  over  this  broad 
frequency  range.  Figure  3.8(b)  shows  the  phase  (ZFn)  over  the  same  range.  As  can  be 
seen  the  agreement  with  the  measured  values  and  the  interpolated  values  is  excellent. 
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Figure  3.6;  Generation  of  stop  band  response  using  pass  band  data. 


(a)  Magnitude  response:  IKnl 
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Figure  3.8:  The  Cauchy  method  applied  to  Device  Characterization. 
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4.0  Effect  of  Noise  in  the  Data  on  the  Cauchy 
Method 


The  Cauchy  method  [7]  has  been  shown  to  provide  accurate  broadband  informa 
tion  from  narrowband  data.  The  method  deals  with  approximating  a  function  by  a 
ratio  of  two  polynomials.  Given  the  value  of  the  function  and  its  derivatives  at  a 
few  points,  the  order  of  the  polynomials  and  their  coefficients  are  evaluated.  Once 
the  coefficients  of  the  two  polynomials  are  known,  they  can  be  used  to  generate  the 
parameter  over  the  entire  band  of  interest.  The  Cauchy  method  was  shown  to  be 
applicable  to  the  cases  where  the  input  data  was  measured  values  of  the  function  and 
not  its  theoretical  values. 

However,  no  measuring  instrument  is  perfect.  Each  measurement  has,  added  to  the 
signal,  an  unwanted  noise  component.  Reference  [7]  does  not  discuss  how  this  noise 
affects  the  results  from  the  Cauchy  method.  The  presence  of  noise  in  the  data  limits 
the  effectiveness  of  the  method.  In  this  chapter  we  try  to  quantify  the  limitations  of 
the  Cauchy  method  when  the  input  data  is  subject  to  contamination  by  noise. 

Throughout  this  work,  we  assume  that  the  noise  is  additive,  stationary,  zero 
mean,  and  Gaussian.  This  assumption  is  approximately  valid  for  most  measuring 
instruments.  Using  this  assumption,  the  probability  density  function  (PDF)  of  the 
parameter,  as  a  function  of  frequency,  is  evaluated.  This  is  compared  to  the  PDF 
approximated  by  a  computer  numerical  simulation. 

To  make  the  problem  tractable,  certain  simplifying  assumptions  are  necessary. 
This  includes  assuming  that  the  coefficients  of  the  polynomials  in  the  Cauchy  method 
are  independent  random  variables.  As  the  theory  will  show,  this  assumption  is  not 
strictly  true.  However,  the  error  introduced  in  the  PDF  due  to  this  assumption  is 
minimal. 

Another  assumption  made  is  that  the  noise  affects  each  measurement  indepen¬ 
dently.  The  noise  is  also  assumed  to  affect  each  measurement,  on  average,  equally. 
This  means,  that  the  average  power  in  the  noise  in  each  measuremeiit  is  assumed 
constant  over  repeated  measurements.  This  assumption  too  is  approximately  valid 
for  most  measuring  systems. 

Using  these  assumptions  we  derive  the  theoretical  PDF  of  the  estimate  of  the 
parameter  as  a  function  of  frequency.  The  theoretical  PDF  was  verified  by  numerical 
simulations. 
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4.1  The  Effect  of  Noise  on  the  Solution  Vector 


As  shown  in  Chapter  2.0,  the  solution  vector  belongs  to  the  invariant  subspace  that 
is  spanned  by  the  right  singular  vector  [V]p+q+2-  This  singular  vector  is  associated 
with  the  smallest  singular  value.  However,  because  of  the  noise  in  the  data,  the  entries 
of  matrix  C  are  perturbed  from  their  true  values.  Hence,  the  solution  vector  is  also 
perturbed.  We  need  to  quantify  the  perturbation  of  this  subspace. 

Notation  In  this  work,  a  perturbed  parameter  or  matrix  will  be  represented  by  a 
tilde  (  )  above  the  corresponding  unperturbed  parameter  or  matrix. 

4.1.1  Perturbation  of  Invariant  Subspaces 

Let,  7?.  denote  the  set  of  real  numbers,  7^",  the  set  of  real  vectors  of  length  n,  and 
7?,"^P  the  set  of  real  matrices  of  order  n  x  p. 

Consider  an  arbitrary  matrix  A  G  with  P  <  N.  Let  A  =  AA  E,  where  E 

is  the  perturbation  to  the  matrix  A,  and 

f  ai  0  \  1 

W'^]\A][V]=  I  2  t  j 

1  P-\ 

Here  the  figures  below  the  matrix  indicate  the  number  of  columns  in  each  submatrix 
while  the  figures  to  the  side  of  the  matrix  indicate  the  number  columns  in  each 
submatrix.  Also, 

U  =  {ui\U2\Us) 

V  =  (U1IV2) 

Here,  u,  €  R",  U,  €  Us  6  ^''■<1''-'’),  e,  6  V.’’,  and  Vs  S 

CTi  is  the  singular  value  corresponding  to  the  left  singular  vector  ui  and  the  right 
singular  vector  v\.  This  singular  value  can  be  the  one  of  interest  and  not  just  the 
largest  singular  value.  In  the  Cauchy  method  the  singular  value  of  interest  is  the 
smallest  or  the  zero  singular  value.  S2  is  the  diagonal  matrix  with  the  rest  of  the 
singular  values  of  C  as  its  entries.  These  singular  values  can  be  ordered  arbitrarily 
as  long  as  the  columns  of  U  and  V  are  permuted  appropriately  so  as  to  maintain  the 
equality  of  equation  (4.1). 

If, 

/  711  9^2  \ 

[UV[E][V]  =  921  G22  (4.2) 

\  931  G32  / 

where  7u  G  72.,  5^12,5^21,5^31  G  72.^''^  Cj22  G  72.^"^^^“\  and  G32  G  and 

if  ai  is  not  repeated  as  a  singular  value,  then  [8] 

5,  =  r,  +  Vs  {isll  -  h  +  0  (|£:|9  (4.3) 
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where, 


h  =  0'\Q\2  "I"  ^2521 

4.1.2  Perturbation  of  the  Solution  in  the  Cauchy  Method 

In  all  measurements  the  true  value  of  the  measured  parameter  (here  Hi)  is  per¬ 
turbed  by  an  additive  noise  component.  Hence, 

Hi  =  Hi  +  e.- 

where.  Hi  is  the  value  of  Hi  after  it  has  been  perturbed  by  noise  e,. 

In  the  following  discussion  we  are  assuming: 

1)  The  noise  is  only  in  the  measurement  of  the  parameter  [/f(s)],  not  in  the  measure¬ 
ment  of  the  frequency  (s). 

2) Cx  =  0  has  a  solution  which  is  unique  to  within  a  constant. 

a)  .x  =  ui,  with  (Tj  =  0 

b)  cTi  =  0  is  a  simple  singular  value.  This  assumption  is  valid  because  m  the  solution 
procedure  we  made  sure  that  the  rank  of  the  null  space  of  C  is  one. 

3)  H{si)  =  H{si)  +  ei,  are  zero  mean,  Gaussian,  uncorrelated,  and  have  equal 

variances  cr^. 

Using  the  above  notation  for  a  perturbed  matrix  and  equation  (2.15),  we  get 


[c 

]  ? 
\[b 

~  0 

(4.4) 

where. 

[C]  = 

1  Si...sf 

1  S2...S^ 

-Hi 

-H2 

-Hls^.. 
-H2S2 . 

..-H2S^ 

(4.5) 

_  1  Sn 

-Hn 

-Hnsn  . . 

.  -  Hns%  . 

where, 

estimate  of  the  order  of  the  numerator, 
Q=:  estimate  of  the  order  of  the  denominator, 
N=  number  of  sample  points. 

=5^  [C]  =  [C]  +  [E]  (4.6) 


where,  E  is  the  additive  error  to  the  matrix  C  due  to  noise  in  the  data.  Hence, 

[E]  =  [0  I  El]  (4.7) 

where,  [0]  is  a  zero  matrix  of  order  N  x  P  +  I  and 


ei 

eiSi 

eisl . . 

Q  1 

.eisj 

e? 

62^2 

6252  •  • 

(4.8) 

eyv 

e^SN 

ejvs^ .  ■ 

Q 

NxQ+1 
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El  =  — 


'  ei 

0 

0  . 

..  0  1 

■  1 

Si 

4 

s^  1 

... 

0 

62 

0  . 

..  0 

1 

S2 

si 

...  *2 

.  0 

0 

0  . 

. .  _ 

1 

Sn 

sjf 

...  Sjsf  j 

(4.9) 


NxN 


Nx.Q+1 


[umiv] = [umEi]  [  v'; 


Using  the  notation  of  section  3.1, 


t 

v;  1 

n 

f;J 

1 

P  +  Q  + 

u'^Eiv;' 

E1V2 

^  7ii 

1  = 

921 

\  5^31 

P+i 
Q  +  i 


iL  \ 

nr 

J'22 


(4.10) 


(4.11) 


(4.12) 


Since  ui  is  the  solution  of  the  unperturbed  Cauchy  equation 


[cl 


=  0 


and  v'l  is  the  vector  of  the  last  Q  +  l  entries  of  Vi,  v'l  is  the  vector  of  denominator 
coefficients.  Also,  the  singular  value  of  interest  (<ri)  is  zero.  Hence,  in  the  notation 
of  section  3.1, 


h  =  (Tigii ^2921 
=  ^2921 


(4.13) 

(4.14) 


also,  using  equation  (4.3) 

Ui  =  m  +  U2S2 '^21  +  0  (4.15) 

Hence,  912  is  of  no  consequence. 

From  equation  (4.12), 

/  711  \ 

[Uf[Ei][v;]  =  921  (4.16) 

\  5^31  / 

Using  equation  (4.9),  and  the  fact  that  Uj  is  the  vector  of  denominator  coefficients. 


Eiv'l  =  - 


Cl  0  0  ...  0 
0  62  0  ...  0 


0  0  0 


de(si) 

de{s2) 

_  de{sN)  . 

(4.17) 
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where,  =  Zk^ohs-,  is  the  value  of  the  unperturbed  denominator  polyno¬ 

mial  evaluated  at  5,  .  For  convenience  we  define  a  new  vector  e  as 


„  // 

EiVi  =  - 


eide{si) 

62^6(52) 


eArde(sjv) 


(4.18) 


Using  this  equation,  the  fact  that  U  =  [ui\U2\U3],  and  equation  (4.16), 

921  =  f^2  ^ 


(4.19) 


Therefore,  using  equation  (4.15), 

til  =  vi  +  V2^2"UU+0  (4.20) 

Since  the  elements  of  e  are  Gaussian  random  variables,  ui  is,  to  the  first  order  of 
approximation  a  Gaussian  random  vector. 

Now,  using  equation  (4.1)  and  the  fact  that  cti  =  0, 

C  =  U2E2VJ 

=»  y2Sj*Uj  =  (4.21) 

where  is  the  psuedo-inverse  of  C. 

Therefore,  to  the  first  order  of  approximation, 

Uj  =  uj  -|-  C^e  (4.22) 

Using  the  fact  that  is  unperturbed  and  the  noise  is  zero  mean,  the  expectation 
value  of  the  solution  vector  (ui)  is  given  by 

E(ni)  =  ui-l-C+E(e)  (4.23) 

(4.24) 

Here  E  is  the  expectation  operator  and  not  the  error  matrix.  Therefore,  to  the  first 
order  of  approximcLtion,  the  estimator  is  unbiased. 

The  covariance  matrix  of  vi  is  given  by, 

cov{vi)  =  E[(vi  —  wi)(ui  —  vi)^] 

Using  equation  (4.22),  we  have 

E((v.  -  t.i)(6,  -  .,)’■)  =  E[C^-£e’’C+’'| 

=  C+El|e’'|C-l’'  (4.25) 
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Now, 


-  ~T 
e  e 


eide[si) 

62^6(52) 


[eide{si)e2de{s2) .  ■  ■  e^de^Sf^)] 


(4.26) 


eNde{sfj)  J 

Therefore,  the  ij-th  entry  of  this  matrix  is  given  by, 

[e  e^]ij  =  eiejde{si)de{sj) 


(4.27) 


Since  e,-  and  Cj  are  assumed  to  be  zero  mean,  independent,  and  identically  distributed 
with  variance  (7^, 


where, 


E[e,ej  (T  Si^ 


I  0  otherwise 


(4.28) 


E[e  e^]  =  a 


de\si)  0  0 

0  de^{s2)  0 


0 


0  0 


E[{vi-vi){vi-vif]  =  c^C^ 


0 

0 


de^(si)  0  0  . . 

0  de^(s2)  0  •  • 


de^{sN) 

0 
0 


(4.29) 


C'+’’  (4.30) 


0  0  0  . . .  de^{sN)  J 

Letting  =  Qj,  the  autocovariance  of  the  z-th  entry  of  vi  is  given  by 


N 


E|C+ee^C+lii  =  <T’X;c?,*“fe) 

i=l 


(4.31) 


This  is  the  variance  of  the  i-th  entry  in  the  vector  of  coefficients.  Hence,  if  i  <  P  +  1, 
we  are  dealing  with  a  numerator  coefficient,  else  we  are  dealing  with  a  denominator 
coefficient. 

Since,  we  have  solved  a  matrix  equation  in  which  the  elements  of  the  matrix  are 
Gaussian  random  variables,  each  element  of  the  solution  vector  is  a  Gaussian  random 
variable.  Also,  the  numerator  and  denominators  are  linear  combinations  of  the  co¬ 
efficients.  Hence,  the  numerator  and  denominator  axe  Gaussian  random  variables  as 
functions  of  frequency.  Hence,  to  completely  characterize  the  numerator  and  denom¬ 
inator  random  variables,  we  only  need  their  expectation  values  and  variances. 

To  make  this  problem  of  the  ratio  of  two  Gaussians  solvable,  we  have  to  assume 
that  any  two  coefficients  are  independent  of  each  other.  Hence,  the  cross-covariance 
matrix  of  vi  is  assumed  to  be  diagonal. 
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Now, 


(4.32) 


Therefore, 


A(5)  = 

k=0 


E[i(s)l  =  EENs' 


(4.33) 


(4.34) 


ElB(«)l  =  EE(6t]s‘  (4.35) 

k=0 

However,  since  to  the  first  order  of  approximation  the  coefficients  are  unbiased, 


E[A(5)]  =  aks’' 

k=0 


(4.36) 


E1B(s)1  =  E‘‘‘'‘ 

Jt=0 


(4.37) 


Therefore,  the  estimators  for  the  numerator  and  denominator  as  a  function  of 
frequency  are  unbiased.  However,  as  we  will  see,  since  the  ratio  of  two  variables  is 
not  a  linear  function,  this  does  not  mean  the  final  estimator  is  unbiased. 

To  calculate  the  variances  of  the  numerator  and  denominator  as  a  function  of 
frequency,  ^ 

uar[A(s)]  =  var  |  ^  o.s*  (4.38) 

.fc=o  J 


(4.38) 


Using  the  assumption  that  each  coefficient  is  independent  of  the  others, 

p 

uar[A(s)]  =  ^  var{ai)s^' 

k=0 

Therefore,  from  equation  (4.31), 


(4.39) 


Similarly, 


P+l  N 

uar[A(s)]  =  ^  c?  de^(sj) 

i=i  j=i 


P+Q+2  N 

var[B{s)]  =  (7^  de^(sj) 

t=P+2  j=i 


(4.40) 


(4.41) 


Let  N  =  E[A(s)],  D  =  E[5(s)],  =  uar[A(s)],  and  =  uar[B(s)].  Therefore, 

the  problem  has  reduced  to;  Given  the  means  and  variances  of  two  independent 
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Gaussian  random  variables,  what  is  the  PDF  of  their  ratio?  This  problem  has  been 
solved  in  reference  [9]. 

In  the  notation  of  [9],  if  N  and  D  are  independent  Gaussian  random  variables 
with  means  N  and  D  respectively,  and  variances  and  respectively,  and  if 


R  = 


N 

D 


then,  the  probability  density  function  of  R  is  given  by 


fnir)  = 


'{abf 


1 


TT 


Zerf{Z)exp{Z'^)  +  -7= 


(4.42) 


where. 


PNr  +  a‘^D\ 
y/PP  +  a?  ) 


and  the  error  function  is  defined  as 


Hence,  we  have  the  theoretical  PDF  of  the  ratio  of  two  random  variables.  However, 
this  density  function  is  an  approximation  of  the  true  density  function.  To  obtain  the 
true  density  function  we  would  need  to  take  into  account  the  cross-correlation  between 
the  coefficients.  This  leads  to  a  problem  that  is  highly  difficult  to  solve. 


4.2  Numerical  Examples 


To  test  the  above  theory,  the  Cauchy  method  was  tested  with  a  simple  example. 
As  an  example  the  function  was  chosen  to  be  the  testing  function  was 


H{s) 


T.U 

ELo(<:+1)»‘ 


(4.43) 


This  ratio  of  two  polynomials  was  evaluated  at  a  31  points  in  the  range  s=2.0  and 
5=4.0.  Two  tests  were  performed  on  this  data. 

In  the  first  test,  Gaussian  noise  was  added  to  the  data  directly.  A  numerical 
Gaussian  random  number  generator  was  used.  The  power  in  the  noise  was  chosen 
such  that  the  signal  to  noise  ratio  {SNR)  was  30dB.  This  perturbed  data  was  used  as 
input  to  the  Cauchy  program.  The  resulting  polynomials  were  used  to  evaluate  the 
parameter  at  s=3.0.  This  was  considered  to  be  one  sample  of  the  random  variable 
at  s=3.0.  1001  such  samples  were  taken.  A  PDF  estimator  was  used  to  estimate  the 
PDF  at  s=3.0.  Figure  4.1  shows  the  PDF  found  using  this  method.  This  is  the  plot 
marked  ’Adding  noise  to  the  data’. 

In  the  second  test,  the  original  unperturbed  data  between  s=2.0  and  s=4.0  were 
used  as  input  to  the  Cauchy  program.  The  unperturbed  numerator  and  denomina¬ 
tor  coefficients  were  evaluated.  The  means  of  the  numerator  and  denominator  were 
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evaluated  using  equations  (4.36)  and  (4.37)  respectively.  Also,  the  variances  of  the 
numerator  and  denominator  were  evaluated  using  equations  (4.40)  and  (4.41)  lespec- 
tively.  Using  these  values  of  means  and  variances,  a  Gaussian  random  variable,  with 
the  numerator  mean  and  variance,  was  divided  with  another  Gaussian  random  vari¬ 
able  with  the  denominator  mean  and  variance.  This  was  repeated  1001  times.  The 
Gaussian  random  numbers  were  generated  using  the  same  random  number  generator 

as  in  the  first  test. 

The  1001  samples  got  from  this  test  were  used  as  input  to  the  same  PDF  estimator. 
The  result  from  this  estimation  of  the  PDF  is  shown  in  Figure  4.1.  This  plot  is  labeled 
’Numerical  Simulation’. 

Finally,  these  two  PDFs  are  compared  with  the  theoretical  PDF  in  equation  (4.42). 
The  choices  of  iV,  D,  a^,  and  6^  are  obtained  from  the  theoretical  means  and  variances 

used  in  the  second  test. 

At  s=3.0,  using  the  above  function 

Actual  Value  :0.2126 
Mean  (adding  noise  to  the  data):0.2124 
Mean  (dividing  two  Gaussians  with  the 
theoretical  means  and  variances): 0.1 804 

Figure  4.2  shows  the  same  three  PDFs  for  a  signal  to  noise  ratio  of  40dB.  Here 
the  agreement  is  better  than  in  the  earlier  case.  This  is  to  be  expected  since  the 
assumptions  come  closer  to  be  being  satisfied  as  the  noise  reduces. 

The  work  in  this  chapter  is  also  published  in  [10] 


5.0  Cauchy  Method  with  Magnitude  Only  Data 


In  Chapter  2.0  the  theory  of  the  Cauchy  method  is  developed  assuming  the  knowl¬ 
edge  of  both  the  real  and  imaginary  parts  of  the  measured  or  numerically  generated 
data  i.e.  the  knowledge  of  the  phase  response  is  assumed  to  be  known.  However,  m 
many  important  applicatons,  it  is  difficult  to  accurately  measure  the  phase  response. 
A  classic  example  of  an  important  situtation  where  this  difficulty  arises  is  the  mea¬ 
suring  the  RCS  of  a  target  in  the  X-band.  The  magnitude  response  can  be  measured 
accurately,  however  calibration  for  the  phase  is  difficult,  if  not  impossible. 

In  this  section  we  propose  possible  approaches  to  extending  the  Ca,uchy  inethod 
to  the  case  where  only  the  magnitude  of  the  data  is  known.  Two  entirely  different 
approaches  are  possible.  In  the  first  approach,  the  Cauchy  method  can  be  used  to 
extrapolate  or  interpolate  the  magnitude  only  data.  The  extrapolated  magnitude 
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only  data  can  then  be  used  as  input  to  a  program  to  retrieve  the  phase  from  the 
magnitude. 

The  second  possible  approach  is  the  reverse  of  the  above  process.  The  magnitude 
only  data  can  be  used  as  input  to  a  program  that  retrieves  the  phase  from  the  given 
magnitude.  This  complex  data  can  then  be  used  by  the  Cauchy  method  described  in 
Chapter  2.0  to  extrapolate  with  respect  to  the  variable  of  interest. 

In  case  of  the  second  approach  the  theory  of  the  Cauchy  method  does  not  require 
any  changes.  The  theory  of  retrieving  the  phase  from  the  magntiude  is  presented  later 
in  this  chapter.  In  the  case  of  the  first  approach,  the  theory  of  the  Cauchy  method 
must  be  modified  to  account  for  the  fact  that  the  given  data  has  only  magnitude 
information. 

5.1  Extrapolation  with  magnitude  only  data 

To  apply  the  Cauchy  method  to  magnitude  only  data,  we  begin  by  recognizing 
that  the  magnitude  response  of  a  linear,  time  invariant,  system  is  an  even  function 
of  frequency  i.e.,  the  magnitude  response  at  fo  is  the  same  as  the  response  at  — /q. 
Hence,  when  modeling  the  magnitude  response  using  a  ratio  of  polynomials,  the  odd 
powers  of  frequency  in  equation  (2.1)  do  not  contribute  to  the  extrapolation.  The 
magnitude  response  \H{s)\  as  a  function  of  frequency  can,  therefore,  be  represented 
as 


i/f(s)i  ^ 


B(s) 


(5.1) 


Since  the  magnitude  response  is  real,  the  coefficients  a{k),k  =  0, ...,P/2  and 
b{k),k  =  0,...,(5/2  that  define  the  two  polynomials  are  real.  Hence  all  the  Her- 
mitian  operations  in  equations  (2.13),  (2.16)  and  (2.17)  must  be  replaced  with  the 
transpose  operator.  Other  than  these  modifications,  the  Cauchy  method,  as  described 
in  Chapter  2.0  is  completely  applicable. 


5.2  Phase  retrieval  from  magnitude  only  data 


Reconstruction  of  the  phase  from  magnitude  only  data  is  an  important  problem. 
For  minimum  phase  systems,  the  reconstruction  of  phase  from  magnitude  only  data 
is  relatively  straightforward  as  the  phase  response  is  given  by  the  Hilbert  transform 
of  the  log  of  the  magnitude  response  [11],  [12],  [13].  Given  the  magnitude  response 
|//’(s)|,  the  phase  response  arg[H{s)]  is  obtained  by 


arg[H{s)]  =  -P 

TT 


/: 


°°  In  \H{x)\ 


dx 


X  —  3 


(5.2) 


where  P  indicates  the  prinicpal  value  integral  as  the  integrand  has  a  singularity  and 
is  not  integrable  in  the  usual  sense. 
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This  relationship  between  the  magnitude  and  the  phase  does  not  hold  if  a  system 
is  not  minimum  phase.  A  system  is  minimum  phase  if  all  zeros  of  the  transfer  function 
are  in  the  left  half  plane  (x  =  a  +  ju;,cr  <  0).  Unfortunately,  most  electromagnetic 
systems  are  non-minimum  phase  and  hence  equation  (5.2)  does  not  apply.  However, 
there  is  a  more  general  result,  based  on  causality,  that  can  be  applied  to  the  phase 
retrieval  problem. 

It  must  be  noted  that  the  phase  retrieval  problem  is  inherently  non-unique.  A  pure 
linear  phase  can  be  added  to  the  reconstructed  phase  without  affecting  the  magnitude 
spectrum.  This  arises  from  the  fact  that  a  linear  phase  is  equivalent  to  a  pure  time 
delay  in  the  time  domain,  the  principle  of  causality  implies  that  the  time  response 
has  to  be  zero  for  negative  times  i.e.,  f{t)  =  0,  f  <  0.  All  physical  systems  must  be 
causal.  A  given  magnitude  response  can  represent  a  causal  system  only  if  it  satisfies 
the  Paley-Wiener  criterion 


r!!M,,<oo  (5.3) 

J-OO  1  +  -S^ 

i.e.  the  integral  must  be  bounded.  Note  that  this  is  possible  only  if  |iJ(s)|  is 
non-zero  for  all  frequencies.  Otherwise,  if  \H{s)\  is  zero  over  any  finite  support,  the 
term  ln|i/(s)l  would  be  infinity.  If  H{$)  =  R{s)  +jX{s),  R{s)  and  A’(s)  are  related 
by  the  Hilbert  transform  [14]  i.e. 


Also, 


and 


.  1  f°°  X(w) 

R(s)  -  —  /  ^  'dw 

TV  J-’OO  S  —  W 

(5.4) 

,  .  1  R{w)  j 

X(s)  =  -  /  dm 

^  ^  TV  J-oo  S  -W 

(5.5) 

h(t)  =  -  /  R(s)  cos(3t)d3 

TV  Jo 

(5.6) 

h{t)  = - /  X{s)siTi{st)ds 

TV  Jo 

(5.7) 

1  roo  ^  TOO 

=  -  |B(^)l'  =  -  / 

TT  J-oo  TT  J-oo 

(5.8) 

In  our  case,  the  magnitude  is  given  over  a  limited  frequency  range  and  at  discrete 
frequency  points.  Hence,  the  integrals  in  the  above  equations  must  be  replaced  by 
summations.  Assuming  that  real  part  R{s)  is  even  with  period  27r  from  —n  to  tt,  we 
can  write  the  Fourier  series. 
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0  <  5  <  TT 


(5.9) 


oo 

R{s)  =  ao  +  XI  a„  cos  (ns  +  4'n) 

n=l 

and  by  the  discrete  version  of  equation  (5.5 


OO 

X(s)  =  —  X  sin(ns  +  <An)  0  <  s  <  tt  (5.10) 

n=l 

i.e.  the  real  and  imaginary  parts  of  the  frequency  response  are  defined  by  the  same 
Fourier  coefficients. 

Now, 


|X(s)p  =  i2^(s)  +  X^(s) 


N 

2 

•  N 

Uo  +  X  cos(ns  +  4n) 

71=1 

+ 

X  «n  sin(ns  +  ^„) 

.n=l 

where,  N  is  chosen  as  a  cut  off  to  limit  the  number  of  unknowns  to  a  finite  number 
based  on  the  information  available.  We  can  now  use  an  optimization  routine  to 
minimize  the  error. 


E{s)  =  \X{s)f 


N 


no  +  X  cos(ns  +  (f>n) 


n=l 


N 


n  2  ' 


X  On  sin(ns  +  <f>n) 


.n=l 


(5.12) 


and  the  phase  can  be  retrieved  from 


-  En=i  On  sin(ns  +  ^„)] 
[oo  +  En=l  On  C0s(ns  + 


(5.13) 


6.0  Conclusions 


This  work  has  presented  a  technique  for  the  determination  of  wideband  response 
of  a  system  given  only  narrowband  information.  Based  on  the  Cauchy  method  starts 
with  assuming  that  the  parameter  of  interest,  as  a  function  of  frequency,  can  be 
approximated  by  a  simple  rational  polynomial  function.  This  assumption  is  valid  for 
many  physical  systems.  The  method  evaluates  the  order  of  the  polynomials  and  the 
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coefficients  that  define  them.  Using  this  form  the  parameter  is  evaluated  at  many 
frequency  points.  It  is  shown  that  the  technique  has  applications  to  many  practical 
problems.  In  our  research  the  technique  is  applied  to  the  Method  of  Moments,  optical 
systems,  filter  analysis,  and  device  characterization.  In  all  applications  the  Cauchy 
method  has  shown  to  save  time  and  memory. 

This  work  has  also  begun  to  explore  the  possibility  of  applying  the  Cauchy  method 
to  magnitude  only  data.  Two  possible  approaches  are  conisdered,  both  of  which  draw 
on  the  ability  to  recover  the  phase  from  the  magnitude  only  data.  A  phase  retrieval 
algorithm  has  been  presented. 

It  must  be  pointed  out  that  the  Cauchy  method  is  completely  general  and  can 
be  used  to  extrapolate  or  interpolate  with  respect  to  any  variable  other  than  fre¬ 
quency.  However,  in  many  applications  in  electromagnetics,  frequency  is  the  variable 

of  interest. 


Bibliography 

[1]  K.  Kottapalli,  T.K.  Sarkar  and  Y.  Hua.  Accurate  computation  of  wide-band 
response  of  electromagnetic  systems  utilizing  narrowband  information.  IEEE 
Transactions  on  Microwave  Theory  and  Techniques,  39(4);682-688,  April  1991. 

[2]  R.F.  Harrington.  Field  Computation  by  Moment  Methods.  Robert  E.  Kreiger 
Publishing  Co.,  1982. 

[3]  G.H.  Golub  and  C.F.  Van  Loan.  Matrix  Computations.  Johns  Hopkins  University 
Press,  1989. 

[4]  S.  Van  Huffel.  Analysis  of  the  Total  Least  Squares  problem  and  its  use  in  param¬ 
eter  estimations.  PhD  thesis,  Kattrolicke  Universiterit  Leuven,  1990. 

[5]  S.M.  Rao.  Electromagnetic  scattering  and  radiation  of  arbitrarily  shaped  surfaces 
by  triangular  patch  modeling.  PhD  thesis.  University  of  Mississipi,  1978. 

[6]  P.W.  Barber  and  S.C.  Hill.  Light  Scattering  by  Particles  -  Computational  Meth¬ 
ods.  World  Scientific,  1990. 


34 


[7]  R.S.  Adve  et.al.  Application  of  the  Cauchy  method  for  extrapolat¬ 
ing/interpolating  narrowband  system  responses.  IEEE  Transactions  on  Mi¬ 
crowave  Theory  and  Techniques,  45,  Part  2(5):837-845,  May  1997. 

[8]  G.W.  Stewart  and  J.  Sun.  Matrix  Perturbation  Theory.  Academic  Press  Inc., 
1990. 

[9]  I.  Kanter.  The  ratios  of  functions  of  random  variables.  IEEE  Transactions  on 
Aerospace  and  Electronic  Systems,  13:624-630,  November  1977. 

[10]  R.  S.  Adve  and  T.  K.  Sarkar.  The  effect  of  noise  in  the  data  on  the  Cauchy 
method.  Microwave  and  Optical  Technology  Letters,  7:242-247,  April  1994. 

[11]  G.C.  James.  Phase  retrieval  and  near  field  antenna  measurements.  In  Proceedings 
of  the  URSI  Symposium,  pages  510-512,  1992.  Sydney,  Australia. 

[12]  Y.M.  Bruck  and  L.G.  Sodin.  On  the  ambiguity  of  the  image  reconstruction 
problem.  Optics  Communications,  30{3):304-308,  September  1979. 

[13]  M.H.  Hayes.  The  reconstruction  of  a  multidimensional  sequence  from  the  phase 
or  magnitude  of  its  Fourier  transform.  IEEE  Transactions  on  Acoustics,  Speech 
and  Signal  Processing,  30(2):140-154,  April  1982. 

[14]  A.V.  Oppenheim  and  R.W.  Schafer.  Discrete-Time  Signal  Processing.  Prentice- 
Hall,  1989. 


35 


MISSION 

OF 

AFRL/INFORMATION DIRECTORATE  (IF) 


The  advancement  and  application  of  information  systems  science  and 
technology  for  aerospace  command  and  control  and  its  transition  to  air, 
space,  and  ground  systems  to  meet  customer  needs  in  the  areas  of  Global 
Awareness,  Dynamic  Planning  and  Execution,  and  Global  Information 
Exchange  is  the  focus  of  this  AFRL  organization.  The  directorate’s  areas 
of  investigation  include  a  broad  spectrum  of  information  and  fusion, 
communication,  collaborative  environment  and  modeling  and  simulation, 
defensive  information  warfare,  and  intelligent  information  systems 
technologies. 


